Goal: Mirror the trajectory analysis performed on the scCLEAN condition
1 - perform Cell Rank according to tutorials: https://cellrank.readthedocs.io/en/stable/index.html
2 - How many lineages are identified?
2 - Do those lineages correlate with coronary vs pulmonary artery?
import scvelo as scv
import scanpy as sc
import cellrank as cr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import seaborn as sns
import scFates as scf
scv.settings.verbosity = 3
scv.settings.set_figure_params("scvelo", dpi_save=600)
cr.settings.verbosity = 2
#open the filtered batch corrected anndata matrix
#samples underwent QC individually, and then aggregated using 'pegasus aggregate_matrix' CLI
#See notebook for identification of Force Directed Layout Trajectory graphs
adata = sc.read_h5ad('VSMC/PC_miQC/PC_miQC_subtypes_Control_VSMC_all_patients.h5ad')
adata
AnnData object with n_obs × n_vars = 34890 × 16168
obs: 'n_genes', 'n_counts', 'percent_mito', 'Channel', 'scale', 'louvain_labels', 'Tissue', 'doublet_score', 'pred_dbl', 'splice', 'value'
var: 'featureid', 'n_cells', 'percent_cells', 'robust', 'highly_variable_features', 'mean', 'var', 'hvf_loess', 'hvf_rank'
uns: 'PCs', '_attr2type', 'diffmap_evals', 'genome', 'louvain_labels_colors', 'louvain_resolution', 'modality', 'ncells', 'norm_count', 'pca', 'pca_features', 'stdzn_max_value', 'stdzn_mean', 'stdzn_std'
obsm: 'X_diffmap', 'X_fle', 'X_pca', 'X_pca_harmony', 'X_phi', 'X_umap', 'diffmap_knn_distances', 'diffmap_knn_indices', 'pca_harmony_knn_distances', 'pca_harmony_knn_indices'
varm: 'means', 'partial_sum'
obsp: 'W_diffmap', 'W_pca_harmony'
adata.obs
| n_genes | n_counts | percent_mito | Channel | scale | louvain_labels | Tissue | doublet_score | pred_dbl | splice | value | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| barcodekey | |||||||||||
| P1_Coronary-AAACCCAAGTCTAGAA | 4458 | 22340 | 7.488809 | P1_Coronary | 4.476276 | 1 | Coronary | 0.001757 | False | Yes | 0 |
| P1_Coronary-AAACCCACAGCACCCA | 6294 | 40694 | 3.351845 | P1_Coronary | 2.457365 | 3 | Coronary | 0.001612 | False | Yes | 0 |
| P1_Coronary-AAACCCAGTTCCGCTT | 3739 | 14339 | 3.528837 | P1_Coronary | 6.974473 | 1 | Coronary | 0.001371 | False | Yes | 0 |
| P1_Coronary-AAACGAAAGAAAGCGA | 4753 | 26268 | 4.301812 | P1_Coronary | 3.807058 | 2 | Coronary | 0.002118 | False | Yes | 0 |
| P1_Coronary-AAACGAAAGGTGGGTT | 3288 | 13315 | 3.604957 | P1_Coronary | 7.510327 | 3 | Coronary | 0.001097 | False | Yes | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| P4_Pulmonary-TTTGGTTAGACTTCGT | 6644 | 61325 | 1.125153 | P4_Pulmonary | 1.630683 | 2 | Pulmonary | 0.001096 | False | Yes | 0 |
| P4_Pulmonary-TTTGTTGAGCATCAAA | 6367 | 40162 | 2.445097 | P4_Pulmonary | 2.490040 | 2 | Pulmonary | 0.001804 | False | Yes | 0 |
| P4_Pulmonary-TTTGTTGAGGTCTGGA | 5293 | 35347 | 3.010156 | P4_Pulmonary | 2.829254 | 1 | Pulmonary | 0.004670 | False | Yes | 0 |
| P4_Pulmonary-TTTGTTGGTGACCTGC | 5182 | 30960 | 3.003876 | P4_Pulmonary | 3.229974 | 1 | Pulmonary | 0.000443 | False | Yes | 0 |
| P4_Pulmonary-TTTGTTGTCTACCTTA | 7676 | 61309 | 2.089416 | P4_Pulmonary | 1.631188 | 2 | Pulmonary | 0.000931 | False | Yes | 0 |
34890 rows × 11 columns
#count the proper number of tissues to ensure the number of colors matches
adata.obs['Tissue'].value_counts()
Pulmonary 19297 Coronary 15593 Name: Tissue, dtype: int64
sc.tl.embedding_density(adata, basis='fle', groupby='Tissue')
sc.set_figure_params(figsize=(6,3),frameon=False,dpi_save=300)
sc.pl.embedding_density(
adata,
basis='fle',
key='fle_density_Tissue',
color_map='GnBu',
save='Density_Control_coronary_pulmonary.png')
WARNING: saving figure to file figures/fle_density_Tissue_Density_Control_coronary_pulmonary.png
/Users/jbezney/opt/anaconda3/envs/CELL_RANK/lib/python3.9/site-packages/scanpy/_settings.py:447: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()` IPython.display.set_matplotlib_formats(*ipython_format) /Users/jbezney/opt/anaconda3/envs/CELL_RANK/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:444: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first. pl.colorbar(cax, ax=ax, pad=0.01, fraction=0.08, aspect=30)
adata.uns['Tissue_colors']= np.array(['#BA0900', '#809693'], dtype=object)
#need to use the harmony batch integrated matrix
sc.pp.neighbors(adata, use_rep='X_pca_harmony')
# see Cell Rank tutorial "CellRank beyond RNA velocity"
# we're copying our `.X` matrix into the layers because that's where `scv.tl.moments`
# expects to find counts for imputation
adata.layers["spliced"] = adata.X
adata.layers["unspliced"] = adata.X
#30 and 30
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
computing neighbors
finished (0:00:02) --> added
'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
finished (0:00:57) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
from cellrank.tl.kernels import CytoTRACEKernel
ctk = CytoTRACEKernel(adata)
Computing CytoTRACE score with `16168` genes
Adding `adata.obs['ct_score']`
`adata.obs['ct_pseudotime']`
`adata.obs['ct_num_exp_genes']`
`adata.var['ct_gene_corr']`
`adata.var['ct_correlates']`
`adata.uns['ct_params']`
Finish (0:00:18)
#plotting using scvelo functions
scv.pl.scatter(adata, basis="fle", color="Tissue", legend_loc="right")
#plotting using scvelo functions
scv.pl.scatter(adata, basis="fle", c="louvain_labels", legend_loc="right")
#plotting using scvelo functions
scv.pl.scatter(adata, basis="fle", c="ct_pseudotime", legend_loc="right",cmap='magma',
title='10x-V3: CytoTRACE Pseudotime', dpi=300, save='VSMC/PC_Figures/Control_Pseudotime.png')
saving figure to file VSMC/PC_Figures/Control_Pseudotime.png
ctk.compute_transition_matrix(threshold_scheme="soft", nu=0.3)
Computing transition matrix based on `ct_pseudotime`
100%|███████████████████████████████████████████████████████████████████████████████████████████████| 34890/34890 [00:07<00:00, 4864.46cell/s]
Finish (0:00:17)
<CytoTRACEKernel>
from cellrank.tl.estimators import GPCCA
g_fwd = GPCCA(ctk)
print(g_fwd)
GPCCA[n=34890, kernel=<CytoTRACEKernel[dnorm=False, scheme=soft, b=10.0, nu=0.3]>]
g_fwd.compute_schur(n_components=25, alpha=0.2)
g_fwd.plot_spectrum(real_only=True, title='10x-V3: Eigengap after 2 Eigenvalues', dpi=300,show_all_xticks=False,
save='VSMC/Figures/Schur_spectrum_control.png')
Computing Schur decomposition
Mat Object: 1 MPI processes
type: seqdense
9.9999999999999978e-01 1.5104703446893388e-03 -1.0589381352350191e-03 2.7950273950795579e-02 5.7465896995287716e-03 -2.6981007403120198e-02 -1.3253115379381099e-02 -6.8190724781176497e-04 -2.6218359832608572e-02 -3.3066354637139618e-02 -5.6385359598725114e-03 1.2441219204278382e-02 1.7873711968797199e-02 -2.7108550577416082e-03 6.9743318436813308e-03 -2.4262715375852124e-02 -7.0706238090262205e-02 -2.5899948554389304e-02 -5.7455826297393542e-02 -9.5251526184787544e-04 6.2442456244687344e-03 -1.7232947388803023e-02 8.3140904450904451e-03 1.0946223742786766e-02 -2.1585266902328183e-02 1.2795137726302936e-02 2.4070688306005745e-02 -2.7892400598665958e-02
0.0000000000000000e+00 9.9627336204564065e-01 5.4521970128935350e-03 9.8249454254520645e-04 -2.6367971432676867e-02 -1.0192294136377173e-02 2.4376187853096100e-04 4.5190890936599318e-03 -4.6776283811514763e-03 -1.9275232991329244e-02 6.5829359852476732e-03 -9.0155665327898236e-03 5.7279085374940559e-04 1.8412026023431481e-03 2.3991022827994819e-02 5.6427958936241017e-03 -3.9958310628242268e-03 -6.2546309197645993e-03 -2.5892796843553598e-02 -2.3725333701813622e-02 -2.1431039784099909e-02 -1.5122492873592554e-02 -1.7984925923894343e-04 3.8945076407209691e-03 -1.1261512805184245e-02 1.3967370317709480e-02 4.0139917876200592e-03 -1.1985099383889086e-02
0.0000000000000000e+00 0.0000000000000000e+00 9.8324778875787489e-01 -8.7536686124265232e-04 -7.4549290234413364e-03 -1.8518199376798053e-03 5.8744103599715823e-04 1.9050220793504818e-03 9.2603528888584501e-04 -4.8631991096830130e-03 3.4140223575729383e-03 -6.5934815214324935e-03 -1.0810547788351638e-03 1.4496609039301213e-03 7.3787955764424342e-03 2.7261606675177451e-03 3.2964656703700426e-03 -1.1114521927562715e-04 -5.0941254438301520e-03 -1.1629490234482219e-02 -8.5878277955733981e-03 -4.4138513740445235e-03 -1.0870281142966222e-03 5.0709826635106597e-04 -4.7832693433645527e-03 5.2446725314915623e-03 5.0538103758852350e-04 -1.4454581112029679e-04
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.7281118371582853e-01 -3.4243483946466316e-03 5.6100879434757226e-03 1.1458710087541700e-03 9.9169271132793467e-03 1.4769863460627279e-02 -3.2024693464858776e-02 1.2453537831746632e-02 -7.8050864185200981e-03 -5.6355978550003596e-03 -3.1423887825281051e-03 3.3412439886181693e-02 -2.3104836168758941e-02 -1.4954065828873250e-02 8.0858369544636460e-03 -2.2799534812089885e-02 -1.8261462149155934e-03 -1.7462234997318536e-03 -1.5881001400214204e-02 2.2256300629418285e-02 1.9404651083168405e-02 -6.0907301755591080e-03 -1.2596475454979378e-02 5.0275210019492507e-03 2.7542019367107911e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.6172932137508682e-01 -8.7845109993197369e-03 -6.1183234209011613e-03 2.7683486888075357e-03 -4.8833301761978928e-04 -5.1794753792574078e-03 9.7072073606289423e-03 -5.6372498432818628e-02 -1.1814715755563980e-03 1.0468379782087318e-03 -3.8194964193303651e-03 7.9853912674039528e-03 9.1965814837766113e-03 -2.0766809648706148e-03 8.7569565838912910e-03 -1.0971524871770765e-02 -3.3215076391372389e-02 -3.8349354240162656e-02 2.6253722125541643e-02 4.3057621918894222e-02 3.0019044654752032e-02 9.7806195813078048e-03 1.1206195923855113e-02 -1.3101122027606577e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.4946879335557066e-01 -1.7090980012679654e-04 -1.0245330396447379e-02 -1.8195505975496695e-02 3.6766217943761594e-02 -1.0713035048589626e-03 -1.0706382856721553e-02 2.9890276086685909e-02 3.8056981527571379e-02 -2.2918869448018173e-02 2.1869217212117480e-02 6.1331634767788062e-03 -2.2821924301566973e-02 4.7976218655189085e-03 7.6767847554851387e-03 -2.0634331301406090e-02 2.1577285663641862e-02 1.7888410356856652e-02 -1.5174625671269159e-02 2.8892934827294404e-03 -8.0278714961815878e-03 3.0097383206820566e-02 2.2826542341472869e-03
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.4855679846759411e-01 -1.2862596905909933e-02 -6.2428715908546069e-04 6.2736654924203670e-03 -6.9806715960521997e-03 -4.7159087161985176e-03 -3.9930468894956847e-03 -1.0789901096621040e-02 -5.2467957874303479e-02 -2.0395440728660508e-02 -5.4624990872841356e-03 -7.0795431961561957e-04 -9.6168405801234022e-03 -1.4203691296345612e-02 -7.9488779336637800e-03 4.9613427323571971e-03 4.4431671382889198e-03 1.9096326291047131e-02 -1.8294854954040707e-02 -3.9929091011835002e-03 -8.0514470496648588e-03 4.0451953589315505e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.3259396545442486e-01 5.0664225718467680e-03 -5.3226975367169750e-03 3.1265498405836883e-03 3.0114429570133245e-03 -1.2228612209173356e-03 -5.7086896471761272e-03 3.7555447879614318e-02 -1.3640248583780384e-02 -3.7023152532333309e-03 7.1199075047312117e-03 -1.3515875430544576e-02 -8.3177962296825548e-03 2.0319152375604558e-02 6.1471860783078036e-03 6.7263133994036068e-03 3.5125755728424774e-03 -4.7053449032174420e-03 8.1749609431236245e-03 2.8391250885571758e-02 3.8228862708859911e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.3158694969544364e-01 2.2242113804779071e-02 -1.7787193291740753e-02 -7.2900936945398259e-03 2.0269979633625746e-02 -6.6844659827845331e-03 -9.0457484384050930e-04 -8.5817263543636335e-03 1.3606310431806688e-02 -1.5827126767738290e-02 -2.8581881263804353e-02 1.6613578669786144e-02 -3.1633676743946244e-02 6.6038715541166498e-02 1.6471657918295527e-02 -1.0278054871803461e-02 1.8014980926325704e-02 -1.3560728951056034e-03 4.4823170021986429e-04 3.0886390633533153e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.3091934441598767e-01 1.4223547926352253e-03 -7.4747519181522131e-04 -4.9652263254875553e-02 5.3607535341217386e-02 -2.2779534248062453e-02 -3.3482297464046587e-02 6.4111895820737980e-02 -1.0135346462070391e-02 -1.0861775008494346e-02 3.2400073128244139e-02 -2.2159238280316892e-02 1.1970995344664784e-02 -3.4941868871880916e-02 1.1929969989059696e-02 1.3990175813765771e-02 -4.7129354204493837e-02 9.1803260027871966e-03 5.2856007136041962e-03
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2991898979791554e-01 1.3587882449117615e-03 1.1393210808298901e-02 -3.2747658221302448e-02 5.7516008468623736e-03 -1.8106653729116958e-02 -1.9787651274014297e-02 4.1712972130151911e-02 -9.4350303229051528e-03 -4.8614276959275384e-03 1.1261269689240464e-02 -4.4003551103625679e-03 -3.9008429620648891e-02 2.5707423059467934e-02 1.8671712238869488e-03 4.0434754961587237e-02 -3.2634577319055787e-02 -5.6334776181945850e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2434023999884818e-01 8.4716453730158178e-03 -1.3527031095716900e-02 1.1060045157607910e-02 7.4876374524998580e-03 -1.1551858947957062e-02 8.8025999974750500e-03 -4.2019534810376477e-03 -1.1374750256859918e-02 -2.2527009073594371e-02 -1.9150335689612559e-02 2.6264754049791554e-02 2.9646350339642732e-02 1.1523702629486302e-02 3.4045959764913611e-03 -1.5961648012890284e-02 -1.7314801715350974e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.1010320312826287e-01 -9.3367368129166940e-03 -1.3417063701558190e-03 5.3608466735433249e-02 -6.2745903247317311e-02 -4.2472737402125184e-02 -1.5292080723881928e-02 3.9842058609974070e-04 -2.5047939541090559e-02 -1.0474608058555660e-02 -2.6229788129411971e-02 4.6616788094295536e-03 -1.2815713736169653e-03 -6.0364812103400542e-02 -7.3459676396517961e-03 2.9940608194055512e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.0468483889695872e-01 -2.1542281197337955e-02 2.1786335835019113e-02 2.0392537366547787e-02 -1.5560562422447339e-02 2.4914337227839884e-03 3.7317987348866695e-03 1.8794826479532879e-02 -1.8739265637647724e-02 -5.1018114331322417e-02 -8.5335098852121691e-03 3.3761831483937361e-02 -1.3001025784296417e-03 -3.9055727948641213e-02 1.6685694701737919e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.0373563739105400e-01 -1.0288112854357262e-02 -3.1813185697879033e-02 -9.1366958887808009e-03 -4.0032634892655204e-02 -2.6667754453496686e-02 2.4844528849893847e-02 2.0599116213199879e-02 7.1591493556318183e-03 -8.2513569852320275e-03 -7.3159237525703311e-04 -4.2282238112723546e-02 1.3276309224545799e-02 2.6528015054847393e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.8788208359093845e-01 1.3198667006794949e-03 2.3096728800521176e-03 5.6197966207582233e-03 -6.0187051411862624e-03 3.2586029140755805e-02 -3.2280725891437663e-02 -1.1619993121017169e-02 -2.7839459884860100e-03 -5.0938190149552309e-03 4.3674424707526881e-02 -1.3045673615349831e-02 4.5703692167642126e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.8540604290006852e-01 -5.2134541440334373e-03 5.4269199247401348e-03 7.6339117295173503e-03 -2.6841612115697045e-03 1.8676373216803252e-02 -1.3716188320876102e-02 -7.7179079727982035e-03 -1.1871583928758622e-02 -1.1924843370880008e-02 2.1031912642256220e-02 -3.8951642414570440e-03
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.8281878410465775e-01 5.7374980637376477e-04 1.2291501419757667e-03 -4.5234218247302855e-03 1.3732616690500144e-02 -2.2669615996593024e-02 2.6285271741407368e-02 -1.8175093128818222e-03 1.2755980199534167e-02 1.6738455944954024e-02 1.0072490460355597e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.7476711397366180e-01 -1.3366865117628298e-02 -1.8867374836053711e-03 1.9918392176175529e-02 3.5547592589591425e-03 -1.4017553556386293e-02 2.7515008565390599e-02 -2.2110532635193358e-02 -1.3644807946228297e-02 -9.7217976505553541e-04
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.7333096535089350e-01 -1.0863173778571196e-02 -1.1433120341631534e-02 1.1776996810847274e-02 2.2805372551047851e-02 9.9717020028976290e-03 -3.1145766570008643e-03 1.0393064067801428e-02 1.4518733214855936e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.7042956314455144e-01 -2.2648521206044311e-02 1.7114528820417699e-02 3.0723547412192902e-02 1.7062258250441002e-02 2.2879725100850734e-02 -1.0139789902297101e-02 -1.8119897238074605e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.6016086987201523e-01 1.7539035439345905e-02 7.1252777802771200e-03 6.0251560770741138e-03 -3.3921803409182480e-02 8.5303425146413935e-03 4.1715796998404650e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.5581726137817804e-01 -2.2127911254151956e-03 -2.1804594232151260e-02 -1.4663114190629576e-02 2.3199767463206001e-02 4.1731446761738199e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.5324121383954465e-01 -1.2541240200928335e-02 1.6663315099352170e-02 -1.2834388874028175e-02 -2.8735122943675678e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4482251210635406e-01 3.2124717459147850e-03 -1.2006066397247104e-02 -2.7960591408977150e-03
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4378110563514586e-01 3.2078737172057856e-03 -1.8088168604686499e-03
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4294492729348747e-01 -2.4050511964898632e-04
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 8.4188552489685931e-01
Adding `adata.uns['eigendecomposition_fwd']`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:02)
According to the gap in the real component of the eigenvalues, two mcarostates should be computed
#cluster_key tries to associated the names of macrostates with cluster labels,
#or in this case, tissue origin
g_fwd.compute_macrostates(n_states=2, cluster_key="Tissue")
g_fwd.plot_macrostates(
discrete=True, legend_loc="right", size=100, basis="fle", title='scCLEAN Macrostates: 2',
figsize=[8,8], dpi=400)
Computing `2` macrostates
Adding `.macrostates`
`.macrostates_memberships`
`.coarse_T`
`.coarse_initial_distribution
`.coarse_stationary_distribution`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:00)
#the course grained transition matrix is stable at 1
#that means there is a clear defined seperation of terminal states
g_fwd.plot_coarse_T(figsize=[5,6], title='10x-V3: Course-grained Transition Matrix',
save='VSMC/PC_Figures/Course_Grained_Control.png', dpi=300)
the higher the stationary distance the more terminal end point in pseudotime
Stationary distance is the stationary distirbution of the transition matrix
The smaller the value the shorter the differentiation path
g_fwd.plot_macrostates(discrete=True, same_plot=True, basis="fle", size=50)
g_fwd.plot_macrostates(same_plot=False, basis='fle',
title=['scCLEAN Macrostate 1','scCLEAN Macrostate 2'],
figsize=[8,6], size=100, fontsize=18, dpi=400)
#add a new observale into adata.obs defining the terminal states
g_fwd.set_terminal_states_from_macrostates(names=['Coronary',
'Pulmonary'])
Adding `adata.obs['terminal_states']`
`adata.obs['terminal_states_probs']`
`.terminal_states`
`.terminal_states_probabilities`
`.terminal_states_memberships
Finish`
g_fwd.rename_terminal_states({'Pulmonary':'Lineage_1', 'Coronary':'Lineage_2'})
#this computes the fate probabilities toward those terminal states for every cell
g_fwd.compute_absorption_probabilities( time_to_absorption='all')
g_fwd.plot_absorption_probabilities(same_plot=False, size=50, basis="fle")
Computing absorption probabilities
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 3.12/s] 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 1.49/s]
Adding `adata.obsm['to_terminal_states']`
`adata.obsm['absorption_times_fwd']`
`.absorption_probabilities`
`.absorption_times`
Finish (0:00:01)
Now it makes sense why the coronary macrostate had a stationary distance of 0
There are only a handful of cells differentiating toward this terminal state
99% of cells are transitioning toward macrostate associated with lineage 2
#check the probabilities of each cell moving to each fate
lineage = pd.DataFrame(adata.obsm['to_terminal_states'])
lineage['Lineage_2'] = lineage[0]>0.5
lineage['Lineage_1'] = lineage[1]>0.5
lineage
| 0 | 1 | Lineage_2 | Lineage_1 | |
|---|---|---|---|---|
| 0 | 1.323540e-12 | 0.999660 | False | True |
| 1 | 2.114814e-12 | 0.999660 | False | True |
| 2 | 9.830771e-13 | 0.999660 | False | True |
| 3 | 4.422072e-12 | 0.999660 | False | True |
| 4 | 1.075595e-12 | 0.999660 | False | True |
| ... | ... | ... | ... | ... |
| 34885 | 8.360096e-10 | 0.999660 | False | True |
| 34886 | 5.840236e-07 | 0.999660 | False | True |
| 34887 | 1.398509e-09 | 0.999660 | False | True |
| 34888 | 1.321043e-09 | 0.999660 | False | True |
| 34889 | 1.363352e-06 | 0.999659 | False | True |
34890 rows × 4 columns
#only 50 cells are differentiation toward lineage 1
lineage['Lineage_2'].value_counts()
False 34840 True 50 Name: Lineage_2, dtype: int64
adata.obs['Tissue'].value_counts()
Pulmonary 19297 Coronary 15593 Name: Tissue, dtype: int64
g_fwd.plot_absorption_probabilities(same_plot=False, size=50,
basis="fle", rescale_color=[0,1.001],
title=['10x-V3 Lineage 2: Absorption Probability','10x-V3 Lineage 1: Absorption Probability'],
figsize=[8,6], fontsize=18, dpi=400, save='VSMC/PC_Figures/Control_absorption_probabilities.png')
saving figure to file VSMC/PC_Figures/Control_absorption_probabilities.png
#use any kwargs from https://scvelo.readthedocs.io/scvelo.pl.scatter/#scvelo.pl.scatter
g_fwd.plot_absorption_probabilities(same_plot=False, color='Tissue', mode='time', time_key='ct_pseudotime',
size=25, basis="fle", alpha=0.5, figsize=[8,8], fontsize=18, dpi=400,
xlabel='Pseudotime', ylabel='Probability',xlim=[-0.05,1.05],
title=['10x-V3 Lineage 2: Absorption Probability','10x-V3 Lineage 1: Absorption Probability'],
save='VSMC/PC_Figures/Control_probability_Coronary_Pulmonary.png')
saving figure to file VSMC/PC_Figures/Control_probability_Coronary_Pulmonary.png
Similar to scCLEAN condition, it is identifying that one lineage is capturing features specific to coronary
Except the amount of cells it can identify as probably differentiating toward that lineage is negligible
scv.pl.scatter(adata, basis="fle", c="terminal_states", size=100, legend_loc='on data',
title='10x-V3 Terminal States', dpi=300, save='VSMC/PC_Figures/Control_1lineage_Terminal_States.png')
saving figure to file VSMC/PC_Figures/Control_1lineage_Terminal_States.png
What is the classification potential for these two lineage distinctions?
#read terminal state probabilities into dataframe to determine correlation
terminal_df = pd.DataFrame(adata.obsm['to_terminal_states'])
terminal_df = terminal_df.rename(columns={0: "Lineage 1", 1: "Lineage 2"})
terminal_df['barcode']=adata.obs.index.values
terminal_df = terminal_df.set_index('barcode')
terminal_df['Origin']=adata.obs['Tissue']
terminal_df
| Lineage 1 | Lineage 2 | Origin | |
|---|---|---|---|
| barcode | |||
| P1_Coronary-AAACCCAAGTCTAGAA | 1.323540e-12 | 0.999660 | Coronary |
| P1_Coronary-AAACCCACAGCACCCA | 2.114814e-12 | 0.999660 | Coronary |
| P1_Coronary-AAACCCAGTTCCGCTT | 9.830771e-13 | 0.999660 | Coronary |
| P1_Coronary-AAACGAAAGAAAGCGA | 4.422072e-12 | 0.999660 | Coronary |
| P1_Coronary-AAACGAAAGGTGGGTT | 1.075595e-12 | 0.999660 | Coronary |
| ... | ... | ... | ... |
| P4_Pulmonary-TTTGGTTAGACTTCGT | 8.360096e-10 | 0.999660 | Pulmonary |
| P4_Pulmonary-TTTGTTGAGCATCAAA | 5.840236e-07 | 0.999660 | Pulmonary |
| P4_Pulmonary-TTTGTTGAGGTCTGGA | 1.398509e-09 | 0.999660 | Pulmonary |
| P4_Pulmonary-TTTGTTGGTGACCTGC | 1.321043e-09 | 0.999660 | Pulmonary |
| P4_Pulmonary-TTTGTTGTCTACCTTA | 1.363352e-06 | 0.999659 | Pulmonary |
34890 rows × 3 columns
adata.obs['Pulmonary_lin_abs_prob']=terminal_df['Lineage 1']
adata.obs['Coronary_lin_abs_prob']=terminal_df['Lineage 2']
#change the identity state into a binary filter
terminal_df['P_identity'] = 0
terminal_df.loc[terminal_df['Origin'] == 'Pulmonary', 'P_identity'] = 1
terminal_df['C_identity'] = 0
terminal_df.loc[terminal_df['Origin'] == 'Coronary', 'C_identity'] = 1
#create new dataframe to generate the heatmap
heat = pd.DataFrame()
heat['Lineage 1'] = terminal_df['Lineage 2']
heat['Pulmonary Identity'] = terminal_df['P_identity']
heat['Lineage 2'] = terminal_df['Lineage 1']
heat['Coronary Identity'] = terminal_df['C_identity']
heat = heat.reset_index(drop=True)
heat
| Lineage 1 | Pulmonary Identity | Lineage 2 | Coronary Identity | |
|---|---|---|---|---|
| 0 | 0.999660 | 0 | 1.323540e-12 | 1 |
| 1 | 0.999660 | 0 | 2.114814e-12 | 1 |
| 2 | 0.999660 | 0 | 9.830771e-13 | 1 |
| 3 | 0.999660 | 0 | 4.422072e-12 | 1 |
| 4 | 0.999660 | 0 | 1.075595e-12 | 1 |
| ... | ... | ... | ... | ... |
| 34885 | 0.999660 | 1 | 8.360096e-10 | 0 |
| 34886 | 0.999660 | 1 | 5.840236e-07 | 0 |
| 34887 | 0.999660 | 1 | 1.398509e-09 | 0 |
| 34888 | 0.999660 | 1 | 1.321043e-09 | 0 |
| 34889 | 0.999659 | 1 | 1.363352e-06 | 0 |
34890 rows × 4 columns
#plot the lineage probability compressed along the pseudotime x-axis
fig = plt.figure()
fig.set_figheight(8)
fig.set_figwidth(6)
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax = sc.pl.violin(adata, 'Coronary_lin_abs_prob', groupby='Tissue', stripplot=True, jitter=0.1, show=False, ax=ax)
ax.set_title('Lineage 2')
ax.set_ylabel('Absorption Probability')
ax.set_ylim(0, 1.05)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax2 = sc.pl.violin(adata, 'Pulmonary_lin_abs_prob', groupby='Tissue', stripplot=True, jitter=0.1, show=False, ax=ax2)
ax2.set_title('Lineage 1')
ax2.set_ylabel('Absorption Probability')
ax2.set_ylim(0, 1.05)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)
fig.savefig('VSMC/PC_Figures/Control_Violin_Lineages.png', dpi=400)
lineage2 = heat[['Lineage 1','Pulmonary Identity']]
fig,ax = plt.subplots(figsize=(3,8))
ax = sns.heatmap(lineage2, cmap='Greys', yticklabels = 5000,
)
plt.xticks(rotation = 30)
# set y-axis label
ax.set_ylabel("Cells",fontsize=12)
ax.figure.axes[-1].set_ylabel('Absorption Probability', size=12, rotation=270, labelpad=20)
fig.autofmt_xdate()
plt.title('Lineage 1 = Pulmonary Lineage', y=1.01, fontsize=12)
# save the plot as a file
fig.savefig('VSMC/PC_Figures/Control_Heatmap_Pulmonary.jpg', dpi=300, bbox_inches='tight')
plt.show()
lineage1 = heat[['Lineage 2','Coronary Identity']]
fig,ax = plt.subplots(figsize=(3,8))
ax = sns.heatmap(lineage1, cmap='Reds', yticklabels = 5000,
)
plt.xticks(rotation = 30)
# set y-axis label
ax.set_ylabel("Cells",fontsize=12)
ax.figure.axes[-1].set_ylabel('Absorption Probability', size=12, rotation=270, labelpad=20)
fig.autofmt_xdate()
plt.title('Lineage 2 = Coronary Lineage', y=1.01, fontsize=12)
# save the plot as a file
fig.savefig('VSMC/PC_Figures/Control_Heatmap_Coronary.jpg', dpi=300, bbox_inches='tight')
plt.show()
#this is for the pulmonary lineage alone exclusively
score = np.array(terminal_df['Lineage 2'])
y = np.array(terminal_df['P_identity'])
# false positive rate
fpr = []
# true positive rate
tpr = []
# Iterate thresholds from 0.0, 0.01, ... 1.0
thresholds = np.arange(0.0, 1.01, .01)
# get number of positive and negative examples in the dataset
P = sum(y)
N = len(y) - P
# iterate through all thresholds and determine fraction of true positives
# and false positives found at this threshold
for thresh in thresholds:
FP=0
TP=0
for i in range(len(score)):
if (score[i] > thresh):
if y[i] == 1:
TP = TP + 1
if y[i] == 0:
FP = FP + 1
fpr.append(FP/float(N))
tpr.append(TP/float(P))
import sklearn.metrics as metrics
roc_auc = metrics.auc(fpr, tpr)
fig, ax = plt.subplots()
fig.set_figheight(8)
fig.set_figwidth(6)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.title('ROC: Lineage 1 = Pulmonary Lineage')
plt.plot(fpr, tpr, 'gray', label = 'AUC = %0.2f' % roc_auc)
plt.plot([0, 1], [0, 1],'black',linestyle='dashed', label = 'Random Classifier')
plt.legend(loc = 'lower right')
plt.xlim([-0.02, 1])
plt.ylim([0, 1.05])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.savefig('VSMC/PC_Figures/Control_ROC_pulmonary_lineage.png', dpi=300, bbox_inches='tight')
plt.show()
#this is for the coronary lineage alone exclusively
score = np.array(terminal_df['Lineage 1'])
y = np.array(terminal_df['C_identity'])
# false positive rate
fpr = []
# true positive rate
tpr = []
# Iterate thresholds from 0.0, 0.01, ... 1.0
thresholds = np.arange(0.0, 1.01, .01)
# get number of positive and negative examples in the dataset
P = sum(y)
N = len(y) - P
# iterate through all thresholds and determine fraction of true positives
# and false positives found at this threshold
for thresh in thresholds:
FP=0
TP=0
for i in range(len(score)):
if (score[i] > thresh):
if y[i] == 1:
TP = TP + 1
if y[i] == 0:
FP = FP + 1
fpr.append(FP/float(N))
tpr.append(TP/float(P))
import sklearn.metrics as metrics
roc_auc = metrics.auc(fpr, tpr)
fig, ax = plt.subplots()
fig.set_figheight(8)
fig.set_figwidth(6)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.title('ROC: Lineage 2 = Coronary Lineage')
plt.plot(fpr, tpr, 'red', label = 'AUC = %0.2f' % roc_auc)
plt.plot([0, 1], [0, 1],'black',linestyle='dashed', label = 'Random Classifier')
plt.legend(loc = 'lower right')
plt.xlim([-0.02, 1])
plt.ylim([0, 1.05])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.savefig('VSMC/PC_Figures/Control_ROC_coronary_lineage.png', dpi=300, bbox_inches='tight')
plt.show()
Clearly, there is no ability to classify lineage association with tissue identity
This was predicted because Lineage 1 had stationary distance along the Markov transition matrix of 0
Furthermore, less than 1% of cells had an absorption probability associated with lineage 1
Consequently, macrostates and lineages will be re-calculated
g_fwd.compute_macrostates(n_states=1, cluster_key="Tissue")
g_fwd.plot_macrostates(
discrete=True, legend_loc="right", size=100, basis="fle", title='scCLEAN Macrostates: 1',
figsize=[8,8], dpi=400)
For 1 macrostate, stationary distribution is computed
Computing eigendecomposition of the transition matrix
Adding `adata.uns['eigendecomposition_fwd']`
`.eigendecomposition`
Finish (0:00:02)
Adding `.macrostates`
`.macrostates_memberships`
`.coarse_T`
`.coarse_initial_distribution
`.coarse_stationary_distribution`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:02)
g_fwd.set_terminal_states_from_macrostates(names=['Pulmonary'])
g_fwd.rename_terminal_states({'Pulmonary':'Lineage_1'})
Adding `adata.obs['terminal_states']`
`adata.obs['terminal_states_probs']`
`.terminal_states`
`.terminal_states_probabilities`
`.terminal_states_memberships
Finish`
g_fwd.compute_absorption_probabilities( time_to_absorption='all')
g_fwd.plot_absorption_probabilities(same_plot=False, size=50, basis="fle")
Computing absorption probabilities WARNING: There is only `1` terminal state, all cells will have probability `1` of going there
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2.44/s] 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.00/s]
Adding `adata.obsm['to_terminal_states']`
`adata.obsm['absorption_times_fwd']`
`.absorption_probabilities`
`.absorption_times`
Finish (0:00:00)
scv.pl.scatter(adata, basis="fle", c="terminal_states", size=100, legend_loc='on data',
title='10x-V3 Terminal State', dpi=300, save='VSMC/PC_Figures/Control_1lineage_Terminal_States.png')
saving figure to file VSMC/PC_Figures/Control_1lineage_Terminal_States.png